/*******************************************************************************
 * Copyright (c) 2012 joey.enfield.
 * All rights reserved. This program and the accompanying materials
 * are made available under the terms of the GNU Public License v3.0
 * which accompanies this distribution, and is available at
 * http://www.gnu.org/licenses/gpl.html
 * 
 * Contributors:
 *     joey.enfield - initial API and implementation
 ******************************************************************************/
package ca.uol.aig.fftpack;

/**
 * @author Baoshe Zhang
 * @author Astronomical Instrument Group of University of Lethbridge.
 */
class RealDoubleFFT_Mixed
{

	/*-------------------------------------------------
	 radf2: Real FFT's forward processing of factor 2
	 -------------------------------------------------*/
	void radf2(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset)
	{
		int i, k, ic;
		double ti2, tr2;
		int iw1;
		iw1 = offset;

		for (k = 0; k < l1; k++)
		{
			ch[2 * k * ido] = cc[k * ido] + cc[(k + l1) * ido];
			ch[(2 * k + 1) * ido + ido - 1] = cc[k * ido] - cc[(k + l1) * ido];
		}
		if (ido < 2)
			return;
		if (ido != 2)
		{
			for (k = 0; k < l1; k++)
			{
				for (i = 2; i < ido; i += 2)
				{
					ic = ido - i;
					tr2 = wtable[i - 2 + iw1] * cc[i - 1 + (k + l1) * ido]
							+ wtable[i - 1 + iw1] * cc[i + (k + l1) * ido];
					ti2 = wtable[i - 2 + iw1] * cc[i + (k + l1) * ido]
							- wtable[i - 1 + iw1] * cc[i - 1 + (k + l1) * ido];
					ch[i + 2 * k * ido] = cc[i + k * ido] + ti2;
					ch[ic + (2 * k + 1) * ido] = ti2 - cc[i + k * ido];
					ch[i - 1 + 2 * k * ido] = cc[i - 1 + k * ido] + tr2;
					ch[ic - 1 + (2 * k + 1) * ido] = cc[i - 1 + k * ido] - tr2;
				}
			}
			if (ido % 2 == 1)
				return;
		}
		for (k = 0; k < l1; k++)
		{
			ch[(2 * k + 1) * ido] = -cc[ido - 1 + (k + l1) * ido];
			ch[ido - 1 + 2 * k * ido] = cc[ido - 1 + k * ido];
		}
	}

	/*-------------------------------------------------
	 radb2: Real FFT's backward processing of factor 2
	 -------------------------------------------------*/
	void radb2(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset)
	{
		int i, k, ic;
		double ti2, tr2;
		int iw1 = offset;

		for (k = 0; k < l1; k++)
		{
			ch[k * ido] = cc[2 * k * ido] + cc[ido - 1 + (2 * k + 1) * ido];
			ch[(k + l1) * ido] = cc[2 * k * ido]
					- cc[ido - 1 + (2 * k + 1) * ido];
		}
		if (ido < 2)
			return;
		if (ido != 2)
		{
			for (k = 0; k < l1; ++k)
			{
				for (i = 2; i < ido; i += 2)
				{
					ic = ido - i;
					ch[i - 1 + k * ido] = cc[i - 1 + 2 * k * ido]
							+ cc[ic - 1 + (2 * k + 1) * ido];
					tr2 = cc[i - 1 + 2 * k * ido]
							- cc[ic - 1 + (2 * k + 1) * ido];
					ch[i + k * ido] = cc[i + 2 * k * ido]
							- cc[ic + (2 * k + 1) * ido];
					ti2 = cc[i + (2 * k) * ido] + cc[ic + (2 * k + 1) * ido];
					ch[i - 1 + (k + l1) * ido] = wtable[i - 2 + iw1] * tr2
							- wtable[i - 1 + iw1] * ti2;
					ch[i + (k + l1) * ido] = wtable[i - 2 + iw1] * ti2
							+ wtable[i - 1 + iw1] * tr2;
				}
			}
			if (ido % 2 == 1)
				return;
		}
		for (k = 0; k < l1; k++)
		{
			ch[ido - 1 + k * ido] = 2 * cc[ido - 1 + 2 * k * ido];
			ch[ido - 1 + (k + l1) * ido] = -2 * cc[(2 * k + 1) * ido];
		}
	}

	/*-------------------------------------------------
	 radf3: Real FFT's forward processing of factor 3 
	 -------------------------------------------------*/
	void radf3(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset)
	{
		final double taur = -0.5D;
		final double taui = 0.866025403784439D;
		int i, k, ic;
		double ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3;
		int iw1, iw2;
		iw1 = offset;
		iw2 = iw1 + ido;

		for (k = 0; k < l1; k++)
		{
			cr2 = cc[(k + l1) * ido] + cc[(k + 2 * l1) * ido];
			ch[3 * k * ido] = cc[k * ido] + cr2;
			ch[(3 * k + 2) * ido] = taui
					* (cc[(k + l1 * 2) * ido] - cc[(k + l1) * ido]);
			ch[ido - 1 + (3 * k + 1) * ido] = cc[k * ido] + taur * cr2;
		}
		if (ido == 1)
			return;
		for (k = 0; k < l1; k++)
		{
			for (i = 2; i < ido; i += 2)
			{
				ic = ido - i;
				dr2 = wtable[i - 2 + iw1] * cc[i - 1 + (k + l1) * ido]
						+ wtable[i - 1 + iw1] * cc[i + (k + l1) * ido];
				di2 = wtable[i - 2 + iw1] * cc[i + (k + l1) * ido]
						- wtable[i - 1 + iw1] * cc[i - 1 + (k + l1) * ido];
				dr3 = wtable[i - 2 + iw2] * cc[i - 1 + (k + l1 * 2) * ido]
						+ wtable[i - 1 + iw2] * cc[i + (k + l1 * 2) * ido];
				di3 = wtable[i - 2 + iw2] * cc[i + (k + l1 * 2) * ido]
						- wtable[i - 1 + iw2] * cc[i - 1 + (k + l1 * 2) * ido];
				cr2 = dr2 + dr3;
				ci2 = di2 + di3;
				ch[i - 1 + 3 * k * ido] = cc[i - 1 + k * ido] + cr2;
				ch[i + 3 * k * ido] = cc[i + k * ido] + ci2;
				tr2 = cc[i - 1 + k * ido] + taur * cr2;
				ti2 = cc[i + k * ido] + taur * ci2;
				tr3 = taui * (di2 - di3);
				ti3 = taui * (dr3 - dr2);
				ch[i - 1 + (3 * k + 2) * ido] = tr2 + tr3;
				ch[ic - 1 + (3 * k + 1) * ido] = tr2 - tr3;
				ch[i + (3 * k + 2) * ido] = ti2 + ti3;
				ch[ic + (3 * k + 1) * ido] = ti3 - ti2;
			}
		}
	}

	/*-------------------------------------------------
	 radb3: Real FFT's backward processing of factor 3
	 -------------------------------------------------*/
	void radb3(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset)
	{
		final double taur = -0.5D;
		final double taui = 0.866025403784439D;
		int i, k, ic;
		double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2;
		int iw1, iw2;
		iw1 = offset;
		iw2 = iw1 + ido;

		for (k = 0; k < l1; k++)
		{
			tr2 = 2 * cc[ido - 1 + (3 * k + 1) * ido];
			cr2 = cc[3 * k * ido] + taur * tr2;
			ch[k * ido] = cc[3 * k * ido] + tr2;
			ci3 = 2 * taui * cc[(3 * k + 2) * ido];
			ch[(k + l1) * ido] = cr2 - ci3;
			ch[(k + 2 * l1) * ido] = cr2 + ci3;
		}
		if (ido == 1)
			return;
		for (k = 0; k < l1; k++)
		{
			for (i = 2; i < ido; i += 2)
			{
				ic = ido - i;
				tr2 = cc[i - 1 + (3 * k + 2) * ido]
						+ cc[ic - 1 + (3 * k + 1) * ido];
				cr2 = cc[i - 1 + 3 * k * ido] + taur * tr2;
				ch[i - 1 + k * ido] = cc[i - 1 + 3 * k * ido] + tr2;
				ti2 = cc[i + (3 * k + 2) * ido] - cc[ic + (3 * k + 1) * ido];
				ci2 = cc[i + 3 * k * ido] + taur * ti2;
				ch[i + k * ido] = cc[i + 3 * k * ido] + ti2;
				cr3 = taui
						* (cc[i - 1 + (3 * k + 2) * ido] - cc[ic - 1
								+ (3 * k + 1) * ido]);
				ci3 = taui
						* (cc[i + (3 * k + 2) * ido] + cc[ic + (3 * k + 1)
								* ido]);
				dr2 = cr2 - ci3;
				dr3 = cr2 + ci3;
				di2 = ci2 + cr3;
				di3 = ci2 - cr3;
				ch[i - 1 + (k + l1) * ido] = wtable[i - 2 + iw1] * dr2
						- wtable[i - 1 + iw1] * di2;
				ch[i + (k + l1) * ido] = wtable[i - 2 + iw1] * di2
						+ wtable[i - 1 + iw1] * dr2;
				ch[i - 1 + (k + 2 * l1) * ido] = wtable[i - 2 + iw2] * dr3
						- wtable[i - 1 + iw2] * di3;
				ch[i + (k + 2 * l1) * ido] = wtable[i - 2 + iw2] * di3
						+ wtable[i - 1 + iw2] * dr3;
			}
		}
	}

	/*-------------------------------------------------
	 radf4: Real FFT's forward processing of factor 4
	 -------------------------------------------------*/
	void radf4(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset)
	{
		final double hsqt2 = 0.7071067811865475D;
		int i, k, ic;
		double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
		int iw1, iw2, iw3;
		iw1 = offset;
		iw2 = offset + ido;
		iw3 = iw2 + ido;
		for (k = 0; k < l1; k++)
		{
			tr1 = cc[(k + l1) * ido] + cc[(k + 3 * l1) * ido];
			tr2 = cc[k * ido] + cc[(k + 2 * l1) * ido];
			ch[4 * k * ido] = tr1 + tr2;
			ch[ido - 1 + (4 * k + 3) * ido] = tr2 - tr1;
			ch[ido - 1 + (4 * k + 1) * ido] = cc[k * ido]
					- cc[(k + 2 * l1) * ido];
			ch[(4 * k + 2) * ido] = cc[(k + 3 * l1) * ido] - cc[(k + l1) * ido];
		}
		if (ido < 2)
			return;
		if (ido != 2)
		{
			for (k = 0; k < l1; k++)
			{
				for (i = 2; i < ido; i += 2)
				{
					ic = ido - i;
					cr2 = wtable[i - 2 + iw1] * cc[i - 1 + (k + l1) * ido]
							+ wtable[i - 1 + iw1] * cc[i + (k + l1) * ido];
					ci2 = wtable[i - 2 + iw1] * cc[i + (k + l1) * ido]
							- wtable[i - 1 + iw1] * cc[i - 1 + (k + l1) * ido];
					cr3 = wtable[i - 2 + iw2] * cc[i - 1 + (k + 2 * l1) * ido]
							+ wtable[i - 1 + iw2] * cc[i + (k + 2 * l1) * ido];
					ci3 = wtable[i - 2 + iw2] * cc[i + (k + 2 * l1) * ido]
							- wtable[i - 1 + iw2]
							* cc[i - 1 + (k + 2 * l1) * ido];
					cr4 = wtable[i - 2 + iw3] * cc[i - 1 + (k + 3 * l1) * ido]
							+ wtable[i - 1 + iw3] * cc[i + (k + 3 * l1) * ido];
					ci4 = wtable[i - 2 + iw3] * cc[i + (k + 3 * l1) * ido]
							- wtable[i - 1 + iw3]
							* cc[i - 1 + (k + 3 * l1) * ido];
					tr1 = cr2 + cr4;
					tr4 = cr4 - cr2;
					ti1 = ci2 + ci4;
					ti4 = ci2 - ci4;
					ti2 = cc[i + k * ido] + ci3;
					ti3 = cc[i + k * ido] - ci3;
					tr2 = cc[i - 1 + k * ido] + cr3;
					tr3 = cc[i - 1 + k * ido] - cr3;
					ch[i - 1 + 4 * k * ido] = tr1 + tr2;
					ch[ic - 1 + (4 * k + 3) * ido] = tr2 - tr1;
					ch[i + 4 * k * ido] = ti1 + ti2;
					ch[ic + (4 * k + 3) * ido] = ti1 - ti2;
					ch[i - 1 + (4 * k + 2) * ido] = ti4 + tr3;
					ch[ic - 1 + (4 * k + 1) * ido] = tr3 - ti4;
					ch[i + (4 * k + 2) * ido] = tr4 + ti3;
					ch[ic + (4 * k + 1) * ido] = tr4 - ti3;
				}
			}
			if (ido % 2 == 1)
				return;
		}
		for (k = 0; k < l1; k++)
		{
			ti1 = -hsqt2
					* (cc[ido - 1 + (k + l1) * ido] + cc[ido - 1 + (k + 3 * l1)
							* ido]);
			tr1 = hsqt2
					* (cc[ido - 1 + (k + l1) * ido] - cc[ido - 1 + (k + 3 * l1)
							* ido]);
			ch[ido - 1 + 4 * k * ido] = tr1 + cc[ido - 1 + k * ido];
			ch[ido - 1 + (4 * k + 2) * ido] = cc[ido - 1 + k * ido] - tr1;
			ch[(4 * k + 1) * ido] = ti1 - cc[ido - 1 + (k + 2 * l1) * ido];
			ch[(4 * k + 3) * ido] = ti1 + cc[ido - 1 + (k + 2 * l1) * ido];
		}
	}

	/*-------------------------------------------------
	 radb4: Real FFT's backward processing of factor 4
	 -------------------------------------------------*/
	void radb4(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset)
	{
		final double sqrt2 = 1.414213562373095D;
		int i, k, ic;
		double ci2, ci3, ci4, cr2, cr3, cr4;
		double ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4;
		int iw1, iw2, iw3;
		iw1 = offset;
		iw2 = iw1 + ido;
		iw3 = iw2 + ido;

		for (k = 0; k < l1; k++)
		{
			tr1 = cc[4 * k * ido] - cc[ido - 1 + (4 * k + 3) * ido];
			tr2 = cc[4 * k * ido] + cc[ido - 1 + (4 * k + 3) * ido];
			tr3 = cc[ido - 1 + (4 * k + 1) * ido]
					+ cc[ido - 1 + (4 * k + 1) * ido];
			tr4 = cc[(4 * k + 2) * ido] + cc[(4 * k + 2) * ido];
			ch[k * ido] = tr2 + tr3;
			ch[(k + l1) * ido] = tr1 - tr4;
			ch[(k + 2 * l1) * ido] = tr2 - tr3;
			ch[(k + 3 * l1) * ido] = tr1 + tr4;
		}
		if (ido < 2)
			return;
		if (ido != 2)
		{
			for (k = 0; k < l1; ++k)
			{
				for (i = 2; i < ido; i += 2)
				{
					ic = ido - i;
					ti1 = cc[i + 4 * k * ido] + cc[ic + (4 * k + 3) * ido];
					ti2 = cc[i + 4 * k * ido] - cc[ic + (4 * k + 3) * ido];
					ti3 = cc[i + (4 * k + 2) * ido]
							- cc[ic + (4 * k + 1) * ido];
					tr4 = cc[i + (4 * k + 2) * ido]
							+ cc[ic + (4 * k + 1) * ido];
					tr1 = cc[i - 1 + 4 * k * ido]
							- cc[ic - 1 + (4 * k + 3) * ido];
					tr2 = cc[i - 1 + 4 * k * ido]
							+ cc[ic - 1 + (4 * k + 3) * ido];
					ti4 = cc[i - 1 + (4 * k + 2) * ido]
							- cc[ic - 1 + (4 * k + 1) * ido];
					tr3 = cc[i - 1 + (4 * k + 2) * ido]
							+ cc[ic - 1 + (4 * k + 1) * ido];
					ch[i - 1 + k * ido] = tr2 + tr3;
					cr3 = tr2 - tr3;
					ch[i + k * ido] = ti2 + ti3;
					ci3 = ti2 - ti3;
					cr2 = tr1 - tr4;
					cr4 = tr1 + tr4;
					ci2 = ti1 + ti4;
					ci4 = ti1 - ti4;
					ch[i - 1 + (k + l1) * ido] = wtable[i - 2 + iw1] * cr2
							- wtable[i - 1 + iw1] * ci2;
					ch[i + (k + l1) * ido] = wtable[i - 2 + iw1] * ci2
							+ wtable[i - 1 + iw1] * cr2;
					ch[i - 1 + (k + 2 * l1) * ido] = wtable[i - 2 + iw2] * cr3
							- wtable[i - 1 + iw2] * ci3;
					ch[i + (k + 2 * l1) * ido] = wtable[i - 2 + iw2] * ci3
							+ wtable[i - 1 + iw2] * cr3;
					ch[i - 1 + (k + 3 * l1) * ido] = wtable[i - 2 + iw3] * cr4
							- wtable[i - 1 + iw3] * ci4;
					ch[i + (k + 3 * l1) * ido] = wtable[i - 2 + iw3] * ci4
							+ wtable[i - 1 + iw3] * cr4;
				}
			}
			if (ido % 2 == 1)
				return;
		}
		for (k = 0; k < l1; k++)
		{
			ti1 = cc[(4 * k + 1) * ido] + cc[(4 * k + 3) * ido];
			ti2 = cc[(4 * k + 3) * ido] - cc[(4 * k + 1) * ido];
			tr1 = cc[ido - 1 + 4 * k * ido] - cc[ido - 1 + (4 * k + 2) * ido];
			tr2 = cc[ido - 1 + 4 * k * ido] + cc[ido - 1 + (4 * k + 2) * ido];
			ch[ido - 1 + k * ido] = tr2 + tr2;
			ch[ido - 1 + (k + l1) * ido] = sqrt2 * (tr1 - ti1);
			ch[ido - 1 + (k + 2 * l1) * ido] = ti2 + ti2;
			ch[ido - 1 + (k + 3 * l1) * ido] = -sqrt2 * (tr1 + ti1);
		}
	}

	/*-------------------------------------------------
	 radf5: Real FFT's forward processing of factor 5
	 -------------------------------------------------*/
	void radf5(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset)
	{
		final double tr11 = 0.309016994374947D;
		final double ti11 = 0.951056516295154D;
		final double tr12 = -0.809016994374947D;
		final double ti12 = 0.587785252292473D;
		int i, k, ic;
		double ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5, cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5;
		int iw1, iw2, iw3, iw4;
		iw1 = offset;
		iw2 = iw1 + ido;
		iw3 = iw2 + ido;
		iw4 = iw3 + ido;

		for (k = 0; k < l1; k++)
		{
			cr2 = cc[(k + 4 * l1) * ido] + cc[(k + l1) * ido];
			ci5 = cc[(k + 4 * l1) * ido] - cc[(k + l1) * ido];
			cr3 = cc[(k + 3 * l1) * ido] + cc[(k + 2 * l1) * ido];
			ci4 = cc[(k + 3 * l1) * ido] - cc[(k + 2 * l1) * ido];
			ch[5 * k * ido] = cc[k * ido] + cr2 + cr3;
			ch[ido - 1 + (5 * k + 1) * ido] = cc[k * ido] + tr11 * cr2 + tr12
					* cr3;
			ch[(5 * k + 2) * ido] = ti11 * ci5 + ti12 * ci4;
			ch[ido - 1 + (5 * k + 3) * ido] = cc[k * ido] + tr12 * cr2 + tr11
					* cr3;
			ch[(5 * k + 4) * ido] = ti12 * ci5 - ti11 * ci4;
		}
		if (ido == 1)
			return;
		for (k = 0; k < l1; ++k)
		{
			for (i = 2; i < ido; i += 2)
			{
				ic = ido - i;
				dr2 = wtable[i - 2 + iw1] * cc[i - 1 + (k + l1) * ido]
						+ wtable[i - 1 + iw1] * cc[i + (k + l1) * ido];
				di2 = wtable[i - 2 + iw1] * cc[i + (k + l1) * ido]
						- wtable[i - 1 + iw1] * cc[i - 1 + (k + l1) * ido];
				dr3 = wtable[i - 2 + iw2] * cc[i - 1 + (k + 2 * l1) * ido]
						+ wtable[i - 1 + iw2] * cc[i + (k + 2 * l1) * ido];
				di3 = wtable[i - 2 + iw2] * cc[i + (k + 2 * l1) * ido]
						- wtable[i - 1 + iw2] * cc[i - 1 + (k + 2 * l1) * ido];
				dr4 = wtable[i - 2 + iw3] * cc[i - 1 + (k + 3 * l1) * ido]
						+ wtable[i - 1 + iw3] * cc[i + (k + 3 * l1) * ido];
				di4 = wtable[i - 2 + iw3] * cc[i + (k + 3 * l1) * ido]
						- wtable[i - 1 + iw3] * cc[i - 1 + (k + 3 * l1) * ido];
				dr5 = wtable[i - 2 + iw4] * cc[i - 1 + (k + 4 * l1) * ido]
						+ wtable[i - 1 + iw4] * cc[i + (k + 4 * l1) * ido];
				di5 = wtable[i - 2 + iw4] * cc[i + (k + 4 * l1) * ido]
						- wtable[i - 1 + iw4] * cc[i - 1 + (k + 4 * l1) * ido];
				cr2 = dr2 + dr5;
				ci5 = dr5 - dr2;
				cr5 = di2 - di5;
				ci2 = di2 + di5;
				cr3 = dr3 + dr4;
				ci4 = dr4 - dr3;
				cr4 = di3 - di4;
				ci3 = di3 + di4;
				ch[i - 1 + 5 * k * ido] = cc[i - 1 + k * ido] + cr2 + cr3;
				ch[i + 5 * k * ido] = cc[i + k * ido] + ci2 + ci3;
				tr2 = cc[i - 1 + k * ido] + tr11 * cr2 + tr12 * cr3;
				ti2 = cc[i + k * ido] + tr11 * ci2 + tr12 * ci3;
				tr3 = cc[i - 1 + k * ido] + tr12 * cr2 + tr11 * cr3;
				ti3 = cc[i + k * ido] + tr12 * ci2 + tr11 * ci3;
				tr5 = ti11 * cr5 + ti12 * cr4;
				ti5 = ti11 * ci5 + ti12 * ci4;
				tr4 = ti12 * cr5 - ti11 * cr4;
				ti4 = ti12 * ci5 - ti11 * ci4;
				ch[i - 1 + (5 * k + 2) * ido] = tr2 + tr5;
				ch[ic - 1 + (5 * k + 1) * ido] = tr2 - tr5;
				ch[i + (5 * k + 2) * ido] = ti2 + ti5;
				ch[ic + (5 * k + 1) * ido] = ti5 - ti2;
				ch[i - 1 + (5 * k + 4) * ido] = tr3 + tr4;
				ch[ic - 1 + (5 * k + 3) * ido] = tr3 - tr4;
				ch[i + (5 * k + 4) * ido] = ti3 + ti4;
				ch[ic + (5 * k + 3) * ido] = ti4 - ti3;
			}
		}
	}

	/*-------------------------------------------------
	 radb5: Real FFT's backward processing of factor 5
	 -------------------------------------------------*/
	void radb5(int ido, int l1, final double cc[], double ch[], final double wtable[], int offset)
	{
		final double tr11 = 0.309016994374947D;
		final double ti11 = 0.951056516295154D;
		final double tr12 = -0.809016994374947D;
		final double ti12 = 0.587785252292473D;
		int i, k, ic;
		double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3, ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5;
		int iw1, iw2, iw3, iw4;
		iw1 = offset;
		iw2 = iw1 + ido;
		iw3 = iw2 + ido;
		iw4 = iw3 + ido;

		for (k = 0; k < l1; k++)
		{
			ti5 = 2 * cc[(5 * k + 2) * ido];
			ti4 = 2 * cc[(5 * k + 4) * ido];
			tr2 = 2 * cc[ido - 1 + (5 * k + 1) * ido];
			tr3 = 2 * cc[ido - 1 + (5 * k + 3) * ido];
			ch[k * ido] = cc[5 * k * ido] + tr2 + tr3;
			cr2 = cc[5 * k * ido] + tr11 * tr2 + tr12 * tr3;
			cr3 = cc[5 * k * ido] + tr12 * tr2 + tr11 * tr3;
			ci5 = ti11 * ti5 + ti12 * ti4;
			ci4 = ti12 * ti5 - ti11 * ti4;
			ch[(k + l1) * ido] = cr2 - ci5;
			ch[(k + 2 * l1) * ido] = cr3 - ci4;
			ch[(k + 3 * l1) * ido] = cr3 + ci4;
			ch[(k + 4 * l1) * ido] = cr2 + ci5;
		}
		if (ido == 1)
			return;
		for (k = 0; k < l1; ++k)
		{
			for (i = 2; i < ido; i += 2)
			{
				ic = ido - i;
				ti5 = cc[i + (5 * k + 2) * ido] + cc[ic + (5 * k + 1) * ido];
				ti2 = cc[i + (5 * k + 2) * ido] - cc[ic + (5 * k + 1) * ido];
				ti4 = cc[i + (5 * k + 4) * ido] + cc[ic + (5 * k + 3) * ido];
				ti3 = cc[i + (5 * k + 4) * ido] - cc[ic + (5 * k + 3) * ido];
				tr5 = cc[i - 1 + (5 * k + 2) * ido]
						- cc[ic - 1 + (5 * k + 1) * ido];
				tr2 = cc[i - 1 + (5 * k + 2) * ido]
						+ cc[ic - 1 + (5 * k + 1) * ido];
				tr4 = cc[i - 1 + (5 * k + 4) * ido]
						- cc[ic - 1 + (5 * k + 3) * ido];
				tr3 = cc[i - 1 + (5 * k + 4) * ido]
						+ cc[ic - 1 + (5 * k + 3) * ido];
				ch[i - 1 + k * ido] = cc[i - 1 + 5 * k * ido] + tr2 + tr3;
				ch[i + k * ido] = cc[i + 5 * k * ido] + ti2 + ti3;
				cr2 = cc[i - 1 + 5 * k * ido] + tr11 * tr2 + tr12 * tr3;

				ci2 = cc[i + 5 * k * ido] + tr11 * ti2 + tr12 * ti3;
				cr3 = cc[i - 1 + 5 * k * ido] + tr12 * tr2 + tr11 * tr3;

				ci3 = cc[i + 5 * k * ido] + tr12 * ti2 + tr11 * ti3;
				cr5 = ti11 * tr5 + ti12 * tr4;
				ci5 = ti11 * ti5 + ti12 * ti4;
				cr4 = ti12 * tr5 - ti11 * tr4;
				ci4 = ti12 * ti5 - ti11 * ti4;
				dr3 = cr3 - ci4;
				dr4 = cr3 + ci4;
				di3 = ci3 + cr4;
				di4 = ci3 - cr4;
				dr5 = cr2 + ci5;
				dr2 = cr2 - ci5;
				di5 = ci2 - cr5;
				di2 = ci2 + cr5;
				ch[i - 1 + (k + l1) * ido] = wtable[i - 2 + iw1] * dr2
						- wtable[i - 1 + iw1] * di2;
				ch[i + (k + l1) * ido] = wtable[i - 2 + iw1] * di2
						+ wtable[i - 1 + iw1] * dr2;
				ch[i - 1 + (k + 2 * l1) * ido] = wtable[i - 2 + iw2] * dr3
						- wtable[i - 1 + iw2] * di3;
				ch[i + (k + 2 * l1) * ido] = wtable[i - 2 + iw2] * di3
						+ wtable[i - 1 + iw2] * dr3;
				ch[i - 1 + (k + 3 * l1) * ido] = wtable[i - 2 + iw3] * dr4
						- wtable[i - 1 + iw3] * di4;
				ch[i + (k + 3 * l1) * ido] = wtable[i - 2 + iw3] * di4
						+ wtable[i - 1 + iw3] * dr4;
				ch[i - 1 + (k + 4 * l1) * ido] = wtable[i - 2 + iw4] * dr5
						- wtable[i - 1 + iw4] * di5;
				ch[i + (k + 4 * l1) * ido] = wtable[i - 2 + iw4] * di5
						+ wtable[i - 1 + iw4] * dr5;
			}
		}
	}

	/*---------------------------------------------------------
	 radfg: Real FFT's forward processing of general factor
	 --------------------------------------------------------*/
	void radfg(int ido, int ip, int l1, int idl1, double cc[], double c1[], double c2[], double ch[], double ch2[], final double wtable[], int offset)
	{
		final double twopi = 2.0D * Math.PI; // 6.28318530717959;
		int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is, nbd;
		double dc2, ai1, ai2, ar1, ar2, ds2, dcp, arg, dsp, ar1h, ar2h;
		int iw1 = offset;

		arg = twopi / ip;
		dcp = Math.cos(arg);
		dsp = Math.sin(arg);
		ipph = (ip + 1) / 2;
		nbd = (ido - 1) / 2;
		if (ido != 1)
		{
			for (ik = 0; ik < idl1; ik++)
				ch2[ik] = c2[ik];
			for (j = 1; j < ip; j++)
				for (k = 0; k < l1; k++)
					ch[(k + j * l1) * ido] = c1[(k + j * l1) * ido];
			if (nbd <= l1)
			{
				is = -ido;
				for (j = 1; j < ip; j++)
				{
					is += ido;
					idij = is - 1;
					for (i = 2; i < ido; i += 2)
					{
						idij += 2;
						for (k = 0; k < l1; k++)
						{
							ch[i - 1 + (k + j * l1) * ido] = wtable[idij - 1
									+ iw1]
									* c1[i - 1 + (k + j * l1) * ido]
									+ wtable[idij + iw1]
									* c1[i + (k + j * l1) * ido];
							ch[i + (k + j * l1) * ido] = wtable[idij - 1 + iw1]
									* c1[i + (k + j * l1) * ido]
									- wtable[idij + iw1]
									* c1[i - 1 + (k + j * l1) * ido];
						}
					}
				}
			} else
			{
				is = -ido;
				for (j = 1; j < ip; j++)
				{
					is += ido;
					for (k = 0; k < l1; k++)
					{
						idij = is - 1;
						for (i = 2; i < ido; i += 2)
						{
							idij += 2;
							ch[i - 1 + (k + j * l1) * ido] = wtable[idij - 1
									+ iw1]
									* c1[i - 1 + (k + j * l1) * ido]
									+ wtable[idij + iw1]
									* c1[i + (k + j * l1) * ido];
							ch[i + (k + j * l1) * ido] = wtable[idij - 1 + iw1]
									* c1[i + (k + j * l1) * ido]
									- wtable[idij + iw1]
									* c1[i - 1 + (k + j * l1) * ido];
						}
					}
				}
			}
			if (nbd >= l1)
			{
				for (j = 1; j < ipph; j++)
				{
					jc = ip - j;
					for (k = 0; k < l1; k++)
					{
						for (i = 2; i < ido; i += 2)
						{
							c1[i - 1 + (k + j * l1) * ido] = ch[i - 1
									+ (k + j * l1) * ido]
									+ ch[i - 1 + (k + jc * l1) * ido];
							c1[i - 1 + (k + jc * l1) * ido] = ch[i
									+ (k + j * l1) * ido]
									- ch[i + (k + jc * l1) * ido];
							c1[i + (k + j * l1) * ido] = ch[i + (k + j * l1)
									* ido]
									+ ch[i + (k + jc * l1) * ido];
							c1[i + (k + jc * l1) * ido] = ch[i - 1
									+ (k + jc * l1) * ido]
									- ch[i - 1 + (k + j * l1) * ido];
						}
					}
				}
			} else
			{
				for (j = 1; j < ipph; j++)
				{
					jc = ip - j;
					for (i = 2; i < ido; i += 2)
					{
						for (k = 0; k < l1; k++)
						{
							c1[i - 1 + (k + j * l1) * ido] = ch[i - 1
									+ (k + j * l1) * ido]
									+ ch[i - 1 + (k + jc * l1) * ido];
							c1[i - 1 + (k + jc * l1) * ido] = ch[i
									+ (k + j * l1) * ido]
									- ch[i + (k + jc * l1) * ido];
							c1[i + (k + j * l1) * ido] = ch[i + (k + j * l1)
									* ido]
									+ ch[i + (k + jc * l1) * ido];
							c1[i + (k + jc * l1) * ido] = ch[i - 1
									+ (k + jc * l1) * ido]
									- ch[i - 1 + (k + j * l1) * ido];
						}
					}
				}
			}
		} else
		{
			for (ik = 0; ik < idl1; ik++)
				c2[ik] = ch2[ik];
		}
		for (j = 1; j < ipph; j++)
		{
			jc = ip - j;
			for (k = 0; k < l1; k++)
			{
				c1[(k + j * l1) * ido] = ch[(k + j * l1) * ido]
						+ ch[(k + jc * l1) * ido];
				c1[(k + jc * l1) * ido] = ch[(k + jc * l1) * ido]
						- ch[(k + j * l1) * ido];
			}
		}

		ar1 = 1;
		ai1 = 0;
		for (l = 1; l < ipph; l++)
		{
			lc = ip - l;
			ar1h = dcp * ar1 - dsp * ai1;
			ai1 = dcp * ai1 + dsp * ar1;
			ar1 = ar1h;
			for (ik = 0; ik < idl1; ik++)
			{
				ch2[ik + l * idl1] = c2[ik] + ar1 * c2[ik + idl1];
				ch2[ik + lc * idl1] = ai1 * c2[ik + (ip - 1) * idl1];
			}
			dc2 = ar1;
			ds2 = ai1;
			ar2 = ar1;
			ai2 = ai1;
			for (j = 2; j < ipph; j++)
			{
				jc = ip - j;
				ar2h = dc2 * ar2 - ds2 * ai2;
				ai2 = dc2 * ai2 + ds2 * ar2;
				ar2 = ar2h;
				for (ik = 0; ik < idl1; ik++)
				{
					ch2[ik + l * idl1] += ar2 * c2[ik + j * idl1];
					ch2[ik + lc * idl1] += ai2 * c2[ik + jc * idl1];
				}
			}
		}
		for (j = 1; j < ipph; j++)
			for (ik = 0; ik < idl1; ik++)
				ch2[ik] += c2[ik + j * idl1];

		if (ido >= l1)
		{
			for (k = 0; k < l1; k++)
			{
				for (i = 0; i < ido; i++)
				{
					cc[i + k * ip * ido] = ch[i + k * ido];
				}
			}
		} else
		{
			for (i = 0; i < ido; i++)
			{
				for (k = 0; k < l1; k++)
				{
					cc[i + k * ip * ido] = ch[i + k * ido];
				}
			}
		}
		for (j = 1; j < ipph; j++)
		{
			jc = ip - j;
			j2 = 2 * j;
			for (k = 0; k < l1; k++)
			{
				cc[ido - 1 + (j2 - 1 + k * ip) * ido] = ch[(k + j * l1) * ido];
				cc[(j2 + k * ip) * ido] = ch[(k + jc * l1) * ido];
			}
		}
		if (ido == 1)
			return;
		if (nbd >= l1)
		{
			for (j = 1; j < ipph; j++)
			{
				jc = ip - j;
				j2 = 2 * j;
				for (k = 0; k < l1; k++)
				{
					for (i = 2; i < ido; i += 2)
					{
						ic = ido - i;
						cc[i - 1 + (j2 + k * ip) * ido] = ch[i - 1
								+ (k + j * l1) * ido]
								+ ch[i - 1 + (k + jc * l1) * ido];
						cc[ic - 1 + (j2 - 1 + k * ip) * ido] = ch[i - 1
								+ (k + j * l1) * ido]
								- ch[i - 1 + (k + jc * l1) * ido];
						cc[i + (j2 + k * ip) * ido] = ch[i + (k + j * l1) * ido]
								+ ch[i + (k + jc * l1) * ido];
						cc[ic + (j2 - 1 + k * ip) * ido] = ch[i + (k + jc * l1)
								* ido]
								- ch[i + (k + j * l1) * ido];
					}
				}
			}
		} else
		{
			for (j = 1; j < ipph; j++)
			{
				jc = ip - j;
				j2 = 2 * j;
				for (i = 2; i < ido; i += 2)
				{
					ic = ido - i;
					for (k = 0; k < l1; k++)
					{
						cc[i - 1 + (j2 + k * ip) * ido] = ch[i - 1
								+ (k + j * l1) * ido]
								+ ch[i - 1 + (k + jc * l1) * ido];
						cc[ic - 1 + (j2 - 1 + k * ip) * ido] = ch[i - 1
								+ (k + j * l1) * ido]
								- ch[i - 1 + (k + jc * l1) * ido];
						cc[i + (j2 + k * ip) * ido] = ch[i + (k + j * l1) * ido]
								+ ch[i + (k + jc * l1) * ido];
						cc[ic + (j2 - 1 + k * ip) * ido] = ch[i + (k + jc * l1)
								* ido]
								- ch[i + (k + j * l1) * ido];
					}
				}
			}
		}
	}

	/*---------------------------------------------------------
	 radbg: Real FFT's backward processing of general factor
	 --------------------------------------------------------*/
	void radbg(int ido, int ip, int l1, int idl1, double cc[], double c1[], double c2[], double ch[], double ch2[], final double wtable[], int offset)
	{
		final double twopi = 2.0D * Math.PI; // 6.28318530717959;
		int idij, ipph, i, j, k, l, j2, ic, jc, lc, ik, is;
		double dc2, ai1, ai2, ar1, ar2, ds2;
		int nbd;
		double dcp, arg, dsp, ar1h, ar2h;
		int iw1 = offset;

		arg = twopi / ip;
		dcp = Math.cos(arg);
		dsp = Math.sin(arg);
		nbd = (ido - 1) / 2;
		ipph = (ip + 1) / 2;
		if (ido >= l1)
		{
			for (k = 0; k < l1; k++)
			{
				for (i = 0; i < ido; i++)
				{
					ch[i + k * ido] = cc[i + k * ip * ido];
				}
			}
		} else
		{
			for (i = 0; i < ido; i++)
			{
				for (k = 0; k < l1; k++)
				{
					ch[i + k * ido] = cc[i + k * ip * ido];
				}
			}
		}
		for (j = 1; j < ipph; j++)
		{
			jc = ip - j;
			j2 = 2 * j;
			for (k = 0; k < l1; k++)
			{
				ch[(k + j * l1) * ido] = cc[ido - 1 + (j2 - 1 + k * ip) * ido]
						+ cc[ido - 1 + (j2 - 1 + k * ip) * ido];
				ch[(k + jc * l1) * ido] = cc[(j2 + k * ip) * ido]
						+ cc[(j2 + k * ip) * ido];
			}
		}

		if (ido != 1)
		{
			if (nbd >= l1)
			{
				for (j = 1; j < ipph; j++)
				{
					jc = ip - j;
					for (k = 0; k < l1; k++)
					{
						for (i = 2; i < ido; i += 2)
						{
							ic = ido - i;
							ch[i - 1 + (k + j * l1) * ido] = cc[i - 1
									+ (2 * j + k * ip) * ido]
									+ cc[ic - 1 + (2 * j - 1 + k * ip) * ido];
							ch[i - 1 + (k + jc * l1) * ido] = cc[i - 1
									+ (2 * j + k * ip) * ido]
									- cc[ic - 1 + (2 * j - 1 + k * ip) * ido];
							ch[i + (k + j * l1) * ido] = cc[i
									+ (2 * j + k * ip) * ido]
									- cc[ic + (2 * j - 1 + k * ip) * ido];
							ch[i + (k + jc * l1) * ido] = cc[i
									+ (2 * j + k * ip) * ido]
									+ cc[ic + (2 * j - 1 + k * ip) * ido];
						}
					}
				}
			} else
			{
				for (j = 1; j < ipph; j++)
				{
					jc = ip - j;
					for (i = 2; i < ido; i += 2)
					{
						ic = ido - i;
						for (k = 0; k < l1; k++)
						{
							ch[i - 1 + (k + j * l1) * ido] = cc[i - 1
									+ (2 * j + k * ip) * ido]
									+ cc[ic - 1 + (2 * j - 1 + k * ip) * ido];
							ch[i - 1 + (k + jc * l1) * ido] = cc[i - 1
									+ (2 * j + k * ip) * ido]
									- cc[ic - 1 + (2 * j - 1 + k * ip) * ido];
							ch[i + (k + j * l1) * ido] = cc[i
									+ (2 * j + k * ip) * ido]
									- cc[ic + (2 * j - 1 + k * ip) * ido];
							ch[i + (k + jc * l1) * ido] = cc[i
									+ (2 * j + k * ip) * ido]
									+ cc[ic + (2 * j - 1 + k * ip) * ido];
						}
					}
				}
			}
		}

		ar1 = 1;
		ai1 = 0;
		for (l = 1; l < ipph; l++)
		{
			lc = ip - l;
			ar1h = dcp * ar1 - dsp * ai1;
			ai1 = dcp * ai1 + dsp * ar1;
			ar1 = ar1h;
			for (ik = 0; ik < idl1; ik++)
			{
				c2[ik + l * idl1] = ch2[ik] + ar1 * ch2[ik + idl1];
				c2[ik + lc * idl1] = ai1 * ch2[ik + (ip - 1) * idl1];
			}
			dc2 = ar1;
			ds2 = ai1;
			ar2 = ar1;
			ai2 = ai1;
			for (j = 2; j < ipph; j++)
			{
				jc = ip - j;
				ar2h = dc2 * ar2 - ds2 * ai2;
				ai2 = dc2 * ai2 + ds2 * ar2;
				ar2 = ar2h;
				for (ik = 0; ik < idl1; ik++)
				{
					c2[ik + l * idl1] += ar2 * ch2[ik + j * idl1];
					c2[ik + lc * idl1] += ai2 * ch2[ik + jc * idl1];
				}
			}
		}
		for (j = 1; j < ipph; j++)
		{
			for (ik = 0; ik < idl1; ik++)
			{
				ch2[ik] += ch2[ik + j * idl1];
			}
		}
		for (j = 1; j < ipph; j++)
		{
			jc = ip - j;
			for (k = 0; k < l1; k++)
			{
				ch[(k + j * l1) * ido] = c1[(k + j * l1) * ido]
						- c1[(k + jc * l1) * ido];
				ch[(k + jc * l1) * ido] = c1[(k + j * l1) * ido]
						+ c1[(k + jc * l1) * ido];
			}
		}

		if (ido == 1)
			return;
		if (nbd >= l1)
		{
			for (j = 1; j < ipph; j++)
			{
				jc = ip - j;
				for (k = 0; k < l1; k++)
				{
					for (i = 2; i < ido; i += 2)
					{
						ch[i - 1 + (k + j * l1) * ido] = c1[i - 1
								+ (k + j * l1) * ido]
								- c1[i + (k + jc * l1) * ido];
						ch[i - 1 + (k + jc * l1) * ido] = c1[i - 1
								+ (k + j * l1) * ido]
								+ c1[i + (k + jc * l1) * ido];
						ch[i + (k + j * l1) * ido] = c1[i + (k + j * l1) * ido]
								+ c1[i - 1 + (k + jc * l1) * ido];
						ch[i + (k + jc * l1) * ido] = c1[i + (k + j * l1) * ido]
								- c1[i - 1 + (k + jc * l1) * ido];
					}
				}
			}
		} else
		{
			for (j = 1; j < ipph; j++)
			{
				jc = ip - j;
				for (i = 2; i < ido; i += 2)
				{
					for (k = 0; k < l1; k++)
					{
						ch[i - 1 + (k + j * l1) * ido] = c1[i - 1
								+ (k + j * l1) * ido]
								- c1[i + (k + jc * l1) * ido];
						ch[i - 1 + (k + jc * l1) * ido] = c1[i - 1
								+ (k + j * l1) * ido]
								+ c1[i + (k + jc * l1) * ido];
						ch[i + (k + j * l1) * ido] = c1[i + (k + j * l1) * ido]
								+ c1[i - 1 + (k + jc * l1) * ido];
						ch[i + (k + jc * l1) * ido] = c1[i + (k + j * l1) * ido]
								- c1[i - 1 + (k + jc * l1) * ido];
					}
				}
			}
		}
		for (ik = 0; ik < idl1; ik++)
			c2[ik] = ch2[ik];
		for (j = 1; j < ip; j++)
			for (k = 0; k < l1; k++)
				c1[(k + j * l1) * ido] = ch[(k + j * l1) * ido];
		if (nbd <= l1)
		{
			is = -ido;
			for (j = 1; j < ip; j++)
			{
				is += ido;
				idij = is - 1;
				for (i = 2; i < ido; i += 2)
				{
					idij += 2;
					for (k = 0; k < l1; k++)
					{
						c1[i - 1 + (k + j * l1) * ido] = wtable[idij - 1 + iw1]
								* ch[i - 1 + (k + j * l1) * ido]
								- wtable[idij + iw1]
								* ch[i + (k + j * l1) * ido];
						c1[i + (k + j * l1) * ido] = wtable[idij - 1 + iw1]
								* ch[i + (k + j * l1) * ido]
								+ wtable[idij + iw1]
								* ch[i - 1 + (k + j * l1) * ido];
					}
				}
			}
		} else
		{
			is = -ido;
			for (j = 1; j < ip; j++)
			{
				is += ido;
				for (k = 0; k < l1; k++)
				{
					idij = is - 1;
					for (i = 2; i < ido; i += 2)
					{
						idij += 2;
						c1[i - 1 + (k + j * l1) * ido] = wtable[idij - 1 + iw1]
								* ch[i - 1 + (k + j * l1) * ido]
								- wtable[idij + iw1]
								* ch[i + (k + j * l1) * ido];
						c1[i + (k + j * l1) * ido] = wtable[idij - 1 + iw1]
								* ch[i + (k + j * l1) * ido]
								+ wtable[idij + iw1]
								* ch[i - 1 + (k + j * l1) * ido];
					}
				}
			}
		}
	}

	/*---------------------------------------------------------
	 rfftf1: further processing of Real forward FFT
	 --------------------------------------------------------*/
	void rfftf1(int n, double c[], final double wtable[], int offset)
	{
		int i;
		int k1, l1, l2, na, kh, nf, ip, iw, ido, idl1;

		double[] ch = new double[n];
		System.arraycopy(wtable, offset, ch, 0, n);

		nf = (int) wtable[1 + 2 * n + offset];
		na = 1;
		l2 = n;
		iw = n - 1 + n + offset;
		for (k1 = 1; k1 <= nf; ++k1)
		{
			kh = nf - k1;
			ip = (int) wtable[kh + 2 + 2 * n + offset];
			l1 = l2 / ip;
			ido = n / l2;
			idl1 = ido * l1;
			iw -= (ip - 1) * ido;
			na = 1 - na;
			if (ip == 4)
			{
				if (na == 0)
				{
					radf4(ido, l1, c, ch, wtable, iw);
				} else
				{
					radf4(ido, l1, ch, c, wtable, iw);
				}
			} else if (ip == 2)
			{
				if (na == 0)
				{
					radf2(ido, l1, c, ch, wtable, iw);
				} else
				{
					radf2(ido, l1, ch, c, wtable, iw);
				}
			} else if (ip == 3)
			{
				if (na == 0)
				{
					radf3(ido, l1, c, ch, wtable, iw);
				} else
				{
					radf3(ido, l1, ch, c, wtable, iw);
				}
			} else if (ip == 5)
			{
				if (na == 0)
				{
					radf5(ido, l1, c, ch, wtable, iw);
				} else
				{
					radf5(ido, l1, ch, c, wtable, iw);
				}
			} else
			{
				if (ido == 1)
					na = 1 - na;
				if (na == 0)
				{
					radfg(ido, ip, l1, idl1, c, c, c, ch, ch, wtable, iw);
					na = 1;
				} else
				{
					radfg(ido, ip, l1, idl1, ch, ch, ch, c, c, wtable, iw);
					na = 0;
				}
			}
			l2 = l1;
		}
		if (na == 1)
			return;
		for (i = 0; i < n; i++)
			c[i] = ch[i];
	}

	/*---------------------------------------------------------
	 rfftb1: further processing of Real backward FFT
	 --------------------------------------------------------*/
	void rfftb1(int n, double c[], final double wtable[], int offset)
	{
		int i;
		int k1, l1, l2, na, nf, ip, iw, ido, idl1;

		double[] ch = new double[n];
		System.arraycopy(wtable, offset, ch, 0, n);

		nf = (int) wtable[1 + 2 * n + offset];
		na = 0;
		l1 = 1;
		iw = n + offset;
		for (k1 = 1; k1 <= nf; k1++)
		{
			ip = (int) wtable[k1 + 1 + 2 * n + offset];
			l2 = ip * l1;
			ido = n / l2;
			idl1 = ido * l1;
			if (ip == 4)
			{
				if (na == 0)
				{
					radb4(ido, l1, c, ch, wtable, iw);
				} else
				{
					radb4(ido, l1, ch, c, wtable, iw);
				}
				na = 1 - na;
			} else if (ip == 2)
			{
				if (na == 0)
				{
					radb2(ido, l1, c, ch, wtable, iw);
				} else
				{
					radb2(ido, l1, ch, c, wtable, iw);
				}
				na = 1 - na;
			} else if (ip == 3)
			{
				if (na == 0)
				{
					radb3(ido, l1, c, ch, wtable, iw);
				} else
				{
					radb3(ido, l1, ch, c, wtable, iw);
				}
				na = 1 - na;
			} else if (ip == 5)
			{
				if (na == 0)
				{
					radb5(ido, l1, c, ch, wtable, iw);
				} else
				{
					radb5(ido, l1, ch, c, wtable, iw);
				}
				na = 1 - na;
			} else
			{
				if (na == 0)
				{
					radbg(ido, ip, l1, idl1, c, c, c, ch, ch, wtable, iw);
				} else
				{
					radbg(ido, ip, l1, idl1, ch, ch, ch, c, c, wtable, iw);
				}
				if (ido == 1)
					na = 1 - na;
			}
			l1 = l2;
			iw += (ip - 1) * ido;
		}
		if (na == 0)
			return;
		for (i = 0; i < n; i++)
			c[i] = ch[i];
	}

	/*---------------------------------------------------------
	 rfftf: Real forward FFT
	 --------------------------------------------------------*/
	void rfftf(int n, double r[], double wtable[])
	{
		if (n == 1)
			return;
		rfftf1(n, r, wtable, 0);
	} /* rfftf */

	/*---------------------------------------------------------
	 rfftf: Real backward FFT
	 --------------------------------------------------------*/
	void rfftb(int n, double r[], double wtable[])
	{
		if (n == 1)
			return;
		rfftb1(n, r, wtable, 0);
	} /* rfftb */

	/*---------------------------------------------------------
	 rffti1: further initialization of Real FFT
	 --------------------------------------------------------*/
	void rffti1(int n, double wtable[], int offset)
	{

		final int[] ntryh = new int[]
		{ 4, 2, 3, 5 };
		final double twopi = 2.0D * Math.PI;
		double argh;
		int ntry = 0, i, j;
		double argld;
		int k1, l1, l2, ib;
		double fi;
		int ld, ii, nf, ip, nl, is, nq, nr;
		double arg;
		int ido, ipm;
		int nfm1;

		nl = n;
		nf = 0;
		j = 0;

		factorize_loop: while (true)
		{
			++j;
			if (j <= 4)
				ntry = ntryh[j - 1];
			else
				ntry += 2;
			do
			{
				nq = nl / ntry;
				nr = nl - ntry * nq;
				if (nr != 0)
					continue factorize_loop;
				++nf;
				wtable[nf + 1 + 2 * n + offset] = ntry;

				nl = nq;
				if (ntry == 2 && nf != 1)
				{
					for (i = 2; i <= nf; i++)
					{
						ib = nf - i + 2;
						wtable[ib + 1 + 2 * n + offset] = wtable[ib + 2 * n
								+ offset];
					}
					wtable[2 + 2 * n + offset] = 2;
				}
			} while (nl != 1);
			break factorize_loop;
		}
		wtable[0 + 2 * n + offset] = n;
		wtable[1 + 2 * n + offset] = nf;
		argh = twopi / (n);
		is = 0;
		nfm1 = nf - 1;
		l1 = 1;
		if (nfm1 == 0)
			return;
		for (k1 = 1; k1 <= nfm1; k1++)
		{
			ip = (int) wtable[k1 + 1 + 2 * n + offset];
			ld = 0;
			l2 = l1 * ip;
			ido = n / l2;
			ipm = ip - 1;
			for (j = 1; j <= ipm; ++j)
			{
				ld += l1;
				i = is;
				argld = ld * argh;

				fi = 0;
				for (ii = 3; ii <= ido; ii += 2)
				{
					i += 2;
					fi += 1;
					arg = fi * argld;
					wtable[i - 2 + n + offset] = Math.cos(arg);
					wtable[i - 1 + n + offset] = Math.sin(arg);
				}
				is += ido;
			}
			l1 = l2;
		}
	} /* rffti1 */

	/*---------------------------------------------------------
	 rffti: Initialization of Real FFT
	 --------------------------------------------------------*/
	void rffti(int n, double wtable[]) /* length of wtable = 2*n + 15 */
	{
		if (n == 1)
			return;
		rffti1(n, wtable, 0);
	} /* rffti */

}
